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Sampling method 



Reld of the Invention 



P cing. 



ai d 



Identifying one or a few states with desired properties, and averaging property s 
of a system, can provide solutions to many important problems including 
binding free energies of drug candidates, protein structure predictions, anafysi 
data, determinations of the value at risk of financial portfolios, or option 
algorithms such as Monte Carlo sampling, Metropolis Monte Carlo aamplfig 
simulations are general methods to address complex optimisation problems 
space integrals. The problem with sampling algorithms is that the number of 
Obtain converged results is often large, and sometimes prohibitively large. The 
With sampling algorithms have, in general, a large number of states x 
probability densities that are known at least up to a normalisation constan 
molecular systems, where x may denote the co-ordinates of the atoms, or 
where x may be a vector of risk factors such as interest rates or stock indices. 



(1) 

(ensity that the 



For Such systems, the average (a) of a property Ok{x) can be calculated witi|Eq. (1), if there 
is a continuum of states, i.e., 
(a)=fp(x)Ok(x)te 

Ih Eq. (1), the integral is over the entire state space, p(x) is the probability 
system is in state x, and the index k may be used to distinguish different propfrties. Similar, if 
the states are discrete, sums 

(OkyZpjOt{xj) (2) 

j 

over all states must be calculated where pj is the probability of the system tcjbe in state xj 
Generalisations to other cases are possible. 

An example is the estimation of the average distance between two atoms frf a molecular 
system. The average distance (R) is equal to the integral 

(R^jcixflixXx) 

over ail conformations of the system. In Eq. (3), x Is the vector of all co-ordinate 
f(x) denotes the probability distribution of the system, /.a, p(x) is equal tc the probability 
density that the system has coordinates x t and £(x) is the distance between tl e two atoms of 
interest as a function of the co-ordinates of the system. 
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of the system, 



'27/08/2002 14:31 

J 

WP-185&-EP 



41-71-911-58-67 



BUEGHEL VON REVY WIL 



08 



10 



15 



20 



25 



30 



•2- 



a state space 
ilobai minimum 
iftegral 

(4) 



The identification of some states with desired properties is a special case o 
integral. E.g.. the identification of the state .*■»(•» that corresponds to the unique 
of some function f{x) that has a unique global minimum is equal to solving the 

In Eq. (4), min(/i(x>)) denotes the unique global minimum of the function f(x). 
Prior Art 



For many problems of interest, the Integrals (Eq. 1) or sums (Eq. 2) cai lot be solved 
analytically, and sampling algorithms have to be used to estimate the props ties. Sampling 
algorithms produce a sample of states such that estimates of the integral of Eq. ( 
Eq. (2) are obtained as a weighted sum over the states of the sample, i.e., 

InEq. (5), Ok denotes the estimate of the average (Ot) of the property Ck(x), / J is the weight 

of the state xtj , the sum is over the N$ states produced by the sampling algorrl im, the index t 
distinguishes different states, the index kcan be used to distinguish different proi erties, and the 
index /can be used to distinguish different runs of the sampling algorithm. 

Established sampling methods include systematic search methods and rar lorn sampling 
methods. Systematic methods, evaluate the properties of the system for s set of states 
determined in advance by a deterministic algorithm. Each state produced neprese its a subset of 
all states of the system. The weight pu used in Eq. (5) for state xu is proportion 1 to the size of 

the subset represented by the state . Systematic search is often used t( estimate low 
dimensional integrals, and improvements of the simple method just described t re known and 
applied routinely, e.g., trapezoidal interpolation. However, there are cases whe e prohibitively 



) or the sum of 



(5) 



large sets of states are required to obtain converged estimates from systematic 
for systems with a large number of degrees of freedom. 



Random sampling methods generate random states xij distributed according o a sampling 
distribution function. For each state considered in run /; the weight pu for E \ (5) can be 
determined, e.g., by Maximum Likelihood considerations (Christian Bartels "An ilysing biased 
Monte Carlo and molecular dynamics simulations 0 , Chem. Phys. Lett 331 (2000 446-454; see 
also M. Souaiile and B. Roux, "Extension to the Weighted Histogram An* ysis Method; 
Combining Umbrella Sampling with Free Energy Calculations", CompuL Phys\ Comm. 135, 



earches, e.g., 



27/08/2002 14:31 



41-71-911-58-57 



'BUECfrEL VON REVY WIL' 



S. 



WP-1859-EP 



-3- 



.30 



' (2001) 40-67; Charles J. Geyer, "Estimating normalizing constants and reweiglilng mixtures' in 
Markov chain Monte Carlo", Tech.. Rep. 568r (1994). School of Statistic! University of 
Minnesota; Art B. Owen and Yi Zhou, "Safe and effective importance samplingl Journal of the 
American Statistical Association 95 (2000) 135-; AM. Ferrenberg and F H. Swendsen, 

5 "Optimized Monte Carlo Data Analysis'; Phys. Rev. Lett. 63, (1 989) 1 1 95-). The ireights depend 
on the sampling distribution functions used. Different random sampling methods such as Monte 
Carlo. Metropolis Monte Carlo, or Langevin Molecular Dynamics differ in the wa ' the states are 
generated. In Monte Carlo sampling - for example according to US 6 061 6 |2 - states are 
generated independently of each other Metropolis Monte Carlo generates new states by 

10 perturbing existing states and accepting or rejecting them such that the stat^ are sampled 
according to a given sampling distribution function. 

With the mentioned sampling algorithms, it is in principle possible to estima s properties of 
arbitrary systems. However, In many cases, the size of the sample required to ol tain converged 
15 . estimates is large, and calculation times get prohibitively long. The convergence Jf the sampling 
algorithm depends critically on the sampling distribution function, which is used. 

For many systems and properties of interest, the best choice of the samp ng distribution 
function is the distribution function of tine system p(x) . This choice of the samf ing distribution 

20 function Is often optimal to estimate the equilibrium properties'' of the system. The weighting 
factors pu of Eq. (5) are all equal to iVf 1 where Ni is the number of states o the sample. In 
other words, estimates of the equilibrium properties of the system are equal to he average of 
the properties of the states sampled. There are cases where the distribution junction of the 
system is not the optimal choice. Examples are the estimation of properties do ninated by the 

25 . contribution of a few rare stales (e.g., the estimation of the probability of a mole ular system to 
be in a transition state region, or the value at risk of a financial portfolio), or syste hs that consist 
of several sub-regfons and where the sampling gets trapped in one of the su -regions (e.g., 
simulations of molecular systems, in general, and Identification of binding modes jf drug-protein 
complexes, in particular). 

The terms "Importance Sampling" and "Umbrella Sampling" refer to sampling llgorithms that 
• use a sampling distribution function different from the distribution function of th< system (MP. 
Allen and D.J. Tildesley, "Computer Simulation of Uquids", Clarendon Press, ( ►xford (1987)). 
These methods are based on the facts that the distribution function of the systei i is not always 
as the optimal sampling distribution function, and that properties of the unbiased ystem can be 
estimated from simulations with a sampling distribution different from the distribu on function of 
the unbiased system. Importance sampling is a general concept. The concep does not say 
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which kind of sampling distribution function to use. For complex systems with a 
states, defining appropriate sampling distribution function is very difficult 



rge number of* 



US .6,381,586 describes a method to get rapid convergence of the estimate of the price of an 
option. The method is based on an appropriate selection of the importance samling distribution 
function. The method relies on the fact that states that contribute to th< ■ estimate are 
concentrated in a sub-region of the state space, and that this sub-region can be i lentified before 
starting, the simulations. If these conditions are fulfilled, convergence of the e itimate can be 
increase by identifying the sub-region, and by defining a sampling distributic } function that 
biases sampling into the identified sub-region. This situation is quite particular ;o the problem 
considered, and the same approach, is not applicable, for example, for the dote nination of the 
value of risk of an option portfolio. In many applications it Is not possible to iden fy such a sub- 
region before starting the simulation. Defining a wrong sub-region would cause t e simulation to 
be trapped within said sub-region and slow down convergence of estimates. 

US 6,178,384 describes a sampling method to estimate conformational free energies. The 
method is based upon an appropriate definition of the sampling distribution functions. The 
method assumes that regions that contribute to the estimates are approximatel>|hanmonic, and 
that the centres of these regions are known. Using these assumptions, a series 
simulations is carried out Each of the simulations samples one of the identified 
sampling distribution function equal to the corresponding harmonic potential, 
difficult to apply, if the system cannot be described as a sum of approximately h*monic wells, if 
there is a large number of harmonic wells that contribute to the estimates, or if tl e wells cannot 
be identified. For example, molecular systems in explicit water consist of a very I irge number of 
wells, or the dominant local minima of proteins are difficult to identify without expe imental data. 



est! nates 



)f Monte Carlo 
sgions using a 
Tie method Is 



Defining sampling distribution function to improve the convergence of the 
sofne knowledge of the system. For most systems of interest, this knowledge is 
priori, and must be acquired somehow. One possibility is to carry through pralimlr fery 
example, minimization as is proposed in the patents 6,178,384 B1 and 6,38 
problem is that with increasing complexity of the system (large number of loca Iminima, 
dimensions, regions of state space that contribute to the estimates are unknot 
preliminary studies become difficult or impossible. 



Instead of defining the sampling distribution function before starting the simulatic X ft is in many 
cases easier to identify a desired property of the sampling, and to iteratively ada t the sampling 
distribution function in order to achieve the desired property. This is the essence of what is 



requires 
not available a 
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586 B1. The 
, many 
i) the required 
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mferred to as "adaptive umbrella sampling* in the literature on molecular simul 



Battels, Michael Schaefer and Martin Karplus determination of equilibria i properties of 



biomolecular systems using multidimensional adaptive umbrella sampling", J. C 



tions (Christian 



em. Phys. 111, 



17 (1999) 8048-8087; Bemd A. Berg and Thomas Neutiaus "Multicanonical e semble: A new 
approach to simulate first-order phase transitions", Phys. Rev. Let 68 (1992), M2). It can for 
example be assumed that the property of interest is the probability distribute i of a selected 
degree of freedom (probability distribution of the distance between two at ms ( probability 
distribution of future gains of a financial institution, ...). Estimates of such a dis ribution require 
that all values of the degree of freedom are sampled Thus, it makes sen? e to adapt the 
sampling distribution function such that all values of interest of the degree )f freedom are 
sampled with comparable probabilities. A problem with adaptive umbrella samp ng is that even 
if the desired umbrella potential is determined, it will in general not be possible 
transitions between different Important regions of the system occur and that the fystern does not 
get trapped 

Summary of the Invention 

The present Invention is a method to improve the performance of establ shed sampling 
methods. The new method can be called incremental umbrella sampling. Increi cental umbrella 
sampling defines a method for sampling the state space by iteratively generate 
their weighting factors pu by fitting the sampling distribution function pj(x) of t 
to at least one weighted property of the already sampled states. This means tfcl pj(x) is fitted 



states xr/ and 
e next iteration 



to the product pjOfas) , in which the pu are the weighting factors and o(xt 
respectively a property of the states xu . 



In a first step an initial sampling distribution function pi(x) is selected. If 
starting parameters are selected as wed, for example a starting state. The iteratidh 
is. set equal to a starting value, preferably to 1. The iteration process consists jjof 
third, a fitting and a fourth step which are repeated until at least one criterion 
second step Nj states xjj are generated by a numerical sampling algorithr 
states according to the sampling distribution function pj(x) . The third step detenfines 
factors pt4 for states xij generated so far f which means for states xij with / 
distinguishes different previous simulations, and the index f distinguishes differer 
simulation). The weighting factors pt are determined by using the sam] 
functions p f (x) used so far, respectively with / <s J . The number of states an< 



) is a function 



nfcessary, other 
parameter / 

a second, a 
fulfilled. In the 

that samples 
weighting 
(the index / 
states of each 
distribution 

the number of 



J 



ip ng 
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Weighting factors pj can be incremented with each iteration. For reason! j of insufficient 
memory ft can be necessary to discard some of the states. In order to have a c insistent set of 
weighting factors p 0 , the weighting factors are calculated in each iteration for ill, respectively 

for a set of seleded, states. At least for normalized weighting factors p if the v\ sighting factors 

5 \ of iterations with / <j (earlier iterations) will be altered by the newest iterati to, respectively 

by the fact that the total number of selected states is increased. The iteration arameter j is 

incremented preferably by 1. The fitting step determines a sampling distribution function p/x) 

: by fitting />X*) to pvOipA for states xtj generated so far. The fitting can be c bne at ail or.at 

, selected states x f/ . In the fourth step at least one criterion is tested, for exampl* the number of 

10 :: iteratiofis. In many cases good results can be found if at least three iterat ins are done.' 
According to the result of the test it will be decided whether to continue the to ration with the 
•7 . second or to go to a filth step in order to perform an analysis. The forth step a bid as welf be 
:l executed directly after the third step. In this case the iteration would be continue a by the fitting 
. step arjd then by the second step. 

15 

By fitting pj{x) in the state space it is possible to use all the information of *nd Oips) for 

: the states x v generated so far. The fitting step allows to use different fitting strategies. For 

* * example the fitting can bias the sampling away from areas where intensive sampling has been 
d.o.ns in the preceding iterations, or the sampling can be directed along ileal gradients 
2b ; • respectively towards local minima or maxima of one or several weighted properips. In each of 
. the iterations, the sampling distribution function is fitted in a way to improve the olbrall sampling 
. ; of the state space. For example each fitted sampling distribution function will dire* the sampling 

• - to* regions that contribute significantly to the estimates of interest and that were litis sampled so 
; . fah For a given region of the state space the fitted sampling distribution is differ Ji)t if there is a 

25 r large number of states with small weighting factors or if there are a few sta es with large 

* .weighting factors. In contrast with adaptive umbrella sampling methods and othe [methods that 
are- based on estimating some properties of the sampling distribution function with =q. 5, the two 
situations would result in comparable sampling distribution functions. 

30. ' Adaptive umbrella sampling attempts to come up with a single sampling distribute i function that 
should ensure sampling of all the important states. Incremental umbrella sampfi ;g creates for 
. each 'iteration a sampling distribution function to improve the overall sampling. Sim * incremental 
umbrella sampling avoids re-sampling of regions that were already extensively sampled and 
'. directs sampling to regions that were relatively little sampled, incremental umbrel p sampling is 

3.5 more efficient than adaptive umbrella sampling. 
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The inventive method reduces the number of states that need to be sample : to identify the 
states of interest, or to obtain converged estimates of the state space integral . Compared to 
simulated annealing (S . Kirkpatrick, C. D. Gelatt Jr., M. P. Vecchi. "Optimizatid i by Simulated 
Annealing", Science, 220, 4598, (1983) 671-680) and multicanonical sampling JBemd A Berg 
and Thomas Neuhaus "Multicanonical ensemble; A new approach to simulate f st-order phase 
transitions", Phys. Rev. Let 68 (1992), 9-12) the invention has several a< vantages. The 
probability of the identified optimal state can be estimated. It supports 
optimisations. State space integrals can be solved It reduces the probability th< 
trapped. The invention is general. It can be used with different sampling methofls, in particular 
With Monte Carlo sampling, Metropolis Monte Carlo sampling, or dynamic simul jb'ons. It can be 



multi-objective 
the system is 



combined with the concepts of simulated annealing and multicanonical sampli 



It provides a 



general framework that can be adapted to the system and the observables of infe lest 

A special embodiment uses in at least one iteration a numerical sampling £ gprithm which 
generates correlated states xj/ , wherein the sampling distribution function pj(x) preferably has 
& maximum that biases sampling, into a sub-region of the state space, where! the sampling 
distribution function pj(x) is preferably fitted such that the maximum is in a re: jon where the 
• product pyOfcj) has a maximum, and wherein the sampling preferably starff from a state 
close to the maximum of the sampling distribution function pj(x) . 

The inventive method Is preferably commercialized In the form of a computer scitware product 
on a computer-readable medium in which program instructions are stored, whkfi instructions, 
when read by a computer, enable the computer to apply the inventive method. 

• Brief Description of the Drawings 

> Fig. 1 . is a flew chart of an iterative umbrella sampling simulation. 
: Fig* 2. shows a distribution function f(x) (asterix), an observable o(x) (circles) 

pfy)0{x) (rectangles) of a one-dimensional system in the form of probabilities 

dimension of the system (x). 

Fig. 3. shows distribution functions pi{x) (circles), pi(x) (rectangles) and pfe)/2+pi(x)f2 
(asterixes) of the one-dimensional system (p versus x). 
. Fig. 4. shows the weighting factors p>\t (circles) and the product pijO(xis) (asteflx) versus the 
states xtj . 



tnd a product 
) versus the 
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f?\jg. 5: is an example of a conformation of the two-dimensional HP proteiivfoldinj model. 

Fig. & shows estimates of the probability (p) of states with a given energy (E) for the HP protein- 

felding model. 

Fig. 7. shows estimates of the probability of states with, energy -92, -90, -88, -86, ;84 and -82 as 
afunction of the number of states that were sampled for the HP protein-folding rr pdel. 
Fig. 8. shows energies (E) of states sampled in test run B as a function of the n imber of states 
(ri) that were sampled. 

Fig. 9; shows estimates of the probability of states (In p) with given radius (R) fbrjuns A, B, C, D 
andE.' 

Fig. to. shows energies (E) of states sampled in the simulated annealing run B is a function of 
this number of states (n) that were sampled. 

Fig. 1 1. shows energies (E) of states sampled in the mulb'canonical (Adumb) run £ as a function 
dfthe number of states (n) that were sampled. 

Fig. 12. shows estimates of the probability (p) of states with a given energy (E). Estimates 
obtained from the five muWcanonical (Adumb) runs A, B, C, D, E and the aveifge of the test 
runs from Fig. 8 Ref. 



: 



25 



35 



Detailed Description of the preferred Embodiments 

Exploration of the state space proceeds in a series of steps shown in Fig. 1. In £ first step 1 an 
initial sampling distribution function pi(x) is selected. If necessary other starting i prameters are 
selected as well, for example a starting state. This is often an interactive step. Computational 
tools may be used. The iteration parameter J is set equal to 1. The iterator .consists of a ■ 
second. 2, a third 3 t a fitting F and a fourth step 4 which are repeated until at least ine criterion is. 
fulfilled. In the second step Nj states xu with t=\^.^Nj are generated fe \ a numerical 
sampling algorithm. In the third step 3 the weighting factors p v are deduced fror | the sampling 
distribution functions p,(x) for states ».# generated so far, which means for st ies xis with / 
£/. Tfje third step 3 increments / by 1. The fitting step F determines the sampl tg distribution 
function pj(x) for the next iteration by fitting pj(x) to psOipj) for states *v ge erated so fan 
Irt the fourth step 4 at least one criterion is tested. For example the criterio i can be the 
convergence of the simulations or reaching a given number of iterations. The criterion may 
invojve some interactive steps. According to the result N or Y of the fourth sb p 4 it will be 
decided whether to continue (N) the iteration with the second step 2 or to go (Y) t a fifth step 5 
In order to perform the final analysis. The simulations, respectively the states xu with the 
weighting factors p itt , are analysed to obtain the estimates of interest An anal sis based on 
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Eq. (5) fa convenient, as It provides a general and consistent way to obt|in estimates of 
equilibrium properties. 



The initial conditions should be selected such that the simulation reaches ra; 
This Is necessary for the sample to represent the sampling distribution functio 
Monte Carlo algorithm or a molecular dynamics simulation is used, the s 
.initial sampling distribution function should be selected such that the starting 
probability of being sampled with the selected sampling distribution function, 
algorithm is used that produces sates independently of any starting state, sele 
state is evidently not required. 



[dly equilibrium. 
If a Metropolis 
state and the 

[ate has a high 
a Monte Carlo 

Ion of a starting 



A simulation is carried out with the selected sampling distribution function /?i(x)|The simulation 
must produce a sample consisting of Nj states xj f distributed according t| the sampling 
distribution function pj(x) . 



Different simulation algorithms may be used, for example Monte Carlo, Metropo^ 
Molecular Dynamics algorithms. The choice will depend on the system to be 
selected sampling distribution function. For reference in later sections,- some 
listed here 

a) Monte Carlo Method: States are generated independently 
transforming random numbers generated by a random number generatoi 
states are un-correlated. 

b) Metropolis Monte Cario, Version 1: Atrial state xm Is 
existing state xy^such that the probability of generating from xjj> 
probability of generating xjj from xww . The trial state is accepted with 

probability 1, if p{x**fcp{xj,) 

probability pfajypfat) if p>{xtr*o)<p(xj/) 

If xtnot is accepted, it is taken as the new state xm=xma , otherwise xu 

new state xjt*\*=x& . 

c) Metropolis Monte Cario, Version 2: Using a Monte Carlo 
produces states distributed according to pm, v i4x), a trial state xww 
independently of any other state. The trial state is accepted with 

probability 1, if p,w{xLf p«Z{xU 



Monte Cario, 
idled and the 
>ssibil'rties are 

each other by 
The resulting 

irated from an 
equal to the 

|{6A) 
;6B) 

taken as the 

Algorithm that 
generated 
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probability /"fr^WM [f ,< -^&L 



(7B) 



If xtrtat Is accepted, it is taken as the new state xjs+i=xtrtai , otherwise xp is tsgt en as the new' 
state xjs+t^x/s . 

New sampling methods are still being developed. US 5,740,072 describes a s|mpling method 
that combines the Metropolis Monte Carlo method with stochastic dynamics. 



hd (9) express 



Optimal weighting factors 

Optimal weighting factors can be obtained with Equations (8-11). Equations (8) 
the weighting factors pff In function of the normalisation constants fff e(f • and v se versa. 



Self-consistent estimates of the weighting factors p'f and normalization const ints Jjf v can 
be obtained by applying equations (8, 9) repeatedly until the weighting factors hq/e converged. 
The normalised weighting factors are then given by 



(9) 



' In equations (8) to (10), Nj and Nk are the number of states produced by simu 
respectively, Nsh* is the number of simulations, and the bias functions 
ct{x)xpk{x)lp(x) 

are: proportional to the ratio of the sampling distribution function of the simulation 
distribution function of the system, it helps with the understanding of^hat follows t 
weighting factors pij are proportional to pfeij)f pfaj) % if only a single simulation 



(10) 

stlons J and k, 
11) 

Kvided by the 
note that the 
s analysed. If 



several simulations are analysed, the weighting factors can be understood as tl p ratio of the 
distribution function of the system divided by an average effective sampling distrit jtion function 
that represents the sampling earned out in all of the simulations. Thus, the weighting factors are 
large for states, for which the distribution function of the system is large, or for whicl the average 
effective sampling distribution function is small. A necessary condition for arw method to 
estimate the weighting factors to be useful in the context of the present mention is that 
additional sampling In the region around a selected state must reduce the weighting factor of the 
selected state. 
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For. particular cases, other equations exist to determine the set of wefghtir ) factors. If the 
sampling distribution functions differ only along one or a few selected degrees :of freedom, the 
Weighted Histogram Analysis Method of A.M. Fenenberg and Swenc [en, "Optimized 
Monte Carlo Data Analysis", Phys. Rev. Lett 63, (1989) 1 195- may be used If t e nonnalfsation 
5 . constants of the sampling distribution functions are known and need not to b< i estimated with 
Eq. (9), the formalisms described and referenced in Art B. Owen and Yi 1 tou, "Safe and 
effective importance sampling 0 , Journal of the American Statistical Associations (2000) 135-' 
may be used. 

jo ; One-dimensional system to iltustrate the inventive method 



A fitting procedure is used to define a sampling distribution function that dire Its sampling to 
regions that contribute significantly to the estimates of interest and that were llle sampled so 
far. One of the advantages of the proposed method is that it can be applied to hph-dlmensionai 
systems with sampling distribution functions that vary along many dimensional The inventive 
method is illustrated using a one-dimensional system, e.g., two atoms separate! by a distance 
X* The distribution function of the system of the example is equal to a normal listribution with 

mean 0 and standard deviation 1 ( i.e. f Mpc} ^ } . The observable ass imed to be of 
interest is set equal to a normal distribution with mean 2 and standard d< Ration 1, Le», 



®h$ h 'r * **^ • Further, it is assumed that two simulations have been carritfi out, and that 



TV 



- the sampling distribution function for the third simulation should be determined, U 

Fig. 2 illustrates the distribution function of the system and the observable; the ai 
distribution function /^(x), the circles the observable Oi(x) , and the rectangles sh 
p(x)CH(x) used in Eq. (1). The product p{x)CX{x) is maximal at x=l. Thus, t< 
averagfe (Ol) of the property Q(x), states around x=l are most important Fij 
sampling distribution functions pi(x) as circles, pz(x) as rectangles, andj 
pfe)/2+pi{x)l2 as asterixes, respectively. The sampling distribution function pi(: 
to be equal to 1 for x between -0.5 and 0.5, the sampling distribution function piif 
to be equal to 0.5 for x between 0 and 2. In each of the two simulations, 3) 
sampled. 

. Fig, 4 plots for each state xy the weighting factors pu and the product pt,tCh(\ 
and asterixes, respectively. The sampling distribution functions are rtormaliSi 
normalisation constants need not to be estimated with Eq. (9), and the (i 



show the 
the product 
calculate the 
3 shdws the 
the average 
was chosen 
was chosen 
states were 

p) as circles 

L Thus, the 
(-normalised) 
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weightingi factors are directly given by 



(12) 



T*iey are equal to the distribution function of the system divided by the average jof the sampling 

distribution functions of the two previous simulations. 
5 The estimate of the average (a) is dominated by states for which the proddbt pijpfaj] Is 

large (Eq. 5). From Eq. (12) and Fig. 4, it is seen that the product pu\a(xu J iJiarge for states 
.• with ffa^afa] large and pifa yi+piipSp. small. p(xi,)a{x>4 is larglfor states that 

contribute significantly to the average (a) of the property Q(x) (Eq. (1 and Fig. 2). 

pi(»j)[2+pi(xj)/2 Is small for states that are In regions that were not much saniled so far (Rg. 
10 ; 3). Thus, to sample states that contribute to the estimate of the average {a}|>f the property 

a(x)i and to sample regions that were little sampled so far, sampling sholld be directed 

towards regions around states that correspond to a large value of the product I^|a(x^J . This 
'. can be done by fitting the sampling distribution function p>(*) for the next sitiulatton to the 

product p,j\Ch{x»} evaluated for the states obtained from the first two simulatons. The fit is 
15 done Such that the resulting sampling distribution function is large for at least sonfe states with a 

large value of the product J»Jpfaj\ . 

1 Formulating this more generally: the sampling distribution function p } (x) for the |ext simulation 
is fitted to the product pvOfa) of the states of all previous simulations times .a 

ap function o{x). The fit is done such that the resulting sampling distribution functiolis large for at 
■ least some states with a large product pvCfov), and tends to be small for Stat* with a small 
' product Ths function C{x) Is defined such that it is large for the sta|s of interest, 

' e^g., it is set equal to c(x)=\a(x\. 

23 Different strategies can be used to fit the sampling distribution function for the nlct iteration to 

• the product pvpifa] shown as asterixes in Rg. 4. A simple strategy is to identiMfhe state with 

• the: largest value of the product A'l^M > and to define *• awnpBng distributirli function for 
the next iteration such that it is large at the selected state. In the example of Fig. 4,|he state with 
thetergest value is the 16* state of the second simulation and is located at 1.07, i.fli, an.w~l.07 . 

30 • Thus, a possible sampling distribution function for the next iteration of the exawfole might be 
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uniform sampling around the state at 1.07 in the range between 0.57 and 1.5* An alternative 
Strategy would be to identify states with a small value of the product ptjp{xtj\ I pd to make the 

sampling distribution function for the next iteration small for these states. B.aL the sampling 
distribution function for the next iteration of the example could be set to a unitirm distribution 
between -0.5 and 2 but excluding the region from -1.37 to 0.63 that lies arounl the state at - 
0.-37 that has the smallest value of the product p<.\a{xv Jj . A third strategy wouli he to express 

the logarithm of the sampling distribution function as a linear function of some parameters and to 
. Use a lineaF least square fit to determine the parameters. These are only somejk the possible 
fitting strategies. They are described in a more formal way below. 1 

Fdrmal description of a fitted sampling distribution function f 

The sampling distribution function p/x) for the next iteration is defined baseclon the states 
xw and their weighting factors pu frpm the previous simulations /=l,2^..,y-l-jrhe sampling 
distribution function pAx) can be expressed as a function h(yj,x) of a vector vj If parameters, 
'. and the normalisation constant fj , i.e., I 

* For some of the existing simulations algorithm, e.g., molecular dynamics or Melopolis Monte 
Carlo, the normalisation constant fj needs not to be known, and needs not to b| determined. 

: For the other simulation algorithms, it can be determined by applying the consiaint that the 
; .integratof the sampling distribution function overall states must be normalised to eye. 

* The' parameters of the sampling distribution function are determined by fitting the sampling 
. dtsftibufion function pj{x) to the weighting factors a» from previous simulations M,2^.,/-l) 

* ' and rejl,2,..,M}. The fit maximises an objective function G. The objective fftiction G Is 
• defined as a function of local comparisons of the sampling distribution function i&ilb*) witn tns 
■: product puOfa) of the weighting factors p f times a function 0{x). The objectrvl function G 
'. must be large for sampling distribution functions that are large for at least some itates With a 
' large product pvOfa), and small for the majority of states with a small produft pyCfa). 

Different procedures can be used to do the fitting; examples are given below. j 

Objective functions G are defined as functions of local comparisons. For the pifoose of the 
•! local comparison, states 35,/ from the previous simulations i=l,2 r 1 are associated to the 
.•weighting factors from previous simulations, e.g., xt^^xi^ . The sampling distribinon function 
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;e function G , 
properties of 

iw. 



• fuiyjJSSs) at x>j is compared to the corresponding product pvOfaj) using a It pal comparison 
. function i 

iiv^(j>j(yj,xs)frs t pi;0(xis)). . (14) 

e.g., <^>(v;^]^ l ^0(%)HM(^}-i>>A'^))' ■ and ^ e objective fur jtion is defined 
5 as. a function of the total comparisons! e.g., <^£ M 2/ ^ ■ 

..; THe form of the sampling distribution function, the set of parameters, the object 
V the function and the fitting procedure may be adapted to the system ar 
v . interest, as well as to the available simulation algorithms. Examples are given bel 

to .; Examtiles of sampling distribution functions with harmonic constraints 

; Examples of sampling distribution functions can be found by applying harmonic jconstraints.fo 
•s the .distribution function of the system, 

15 . In Eq. (15), fcawrr is a constant that determines the strength of the harmonic 
the scaling factor to normalise pj(x) , and xjw«r is the centre of the constraint 
■ that cgn be fitted are the centre and the strength of the constraint, vj^kcar^^xj^m . 
; For systems consisting of many local minima separated ..by barriers, apljication of a 
. muiticanonical sampling distribution function can be advantageous, 

In Eq. {16), pmJ{x) denotes the multicanonical distribution function; the remaining parameters 
. are/identical to those of Eq. (15). In particular, the parameters that can be fitted #e the centre 
■■> antfthe' strength of the constraint, vj={kc 0 >a(r,xj^m J tr} . 



:i5> 

istraint, fi is 
, y parameters 



25 'Th> multicanonical distribution function p*<{x) (Bemd A. Berg and Thorris Neuhaus 
' "Mtfjticarjonical ensemble: A new approach to simulate first-order phase transition*, Phys. Rev. 
*Let; 69 (1992); Christian Bartels, Michael Schaefer and Martin Karplus "Del Irmination of 
• equilibrium properties of biomolecutar systems using multidimensional adapl re umbrella 
Sampling", J. Chem. Phys. 111, 17 (1999) 8048-8067) is determined such that dii isrent values 

30 of ln/j(r) are sampled with comparable probability. The multicanonical distributer function can 
be Written as 

/& IC (jc>cp(xy/n. c (ln^)), ( j 7 ) 

. .where /«c denotes an estimate of the probability distribution of Vap{x) . The estin ste fa can 
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be obtained from a histogram of kip(x) that can be obtained from the sampfc !i states jkv and 
..their weighting factors jev using Eq. (5) . 

A sum of harmonic constraints applied to the distribution function of the syst€ h results in Hie 
distribution function 

♦ M*)=^exjp^ ! (18) 

to Eq.{18)i the constants kj^w determine the strengths of the hannonic consi ^ints, fj is the 
sealing factor to normalise pj(x), xjs,con*r are the centres of the constraints, an| Nj^mur is the 
number of hannonic constraints. The parameters that can be fitted are the numfer of harmonic 
constraints* their centres, and their strengths, 



iods may be 
(15) or (f 6) 

of the state 
jpolis Monte 
itrast, if the 
Monte Carto 
J Monte Carlo 



Depending on the sampling distribution function pj(x), particular sampling mi 
. fridrexjr less useful. For example, the sampling distribution functions defined by 
. with a/positive value for the constant constrain sampling to a sub-regl 
; space and are particular useful together with sampling methods such as the 
* Carlo Version 1 or stochastic dynamics that generate correlated states. In 

sampling distribution function of Eq. (18) is used, convergence of a Metropolis 
: Version 1 algorithm may be slow, and Metropolis Monte Carto Version 2 oi 
. algorithms may be preferable. 

*i • 

. Examples of the function cXx) 



Examples of the function 0(x) can be defined such that p{x)o(x) is large for stat |s that should 
. b&sampled in one of the simulations. The states that should be sampled are prln krily those for 
which one of the N P rc P properties 0^a(xlOi(x}... 9 O)rw(x)} (see below and E }. 1) is.largfe 
The>funetion o(x) is set equal to a function of the set 0 of properties, e.g., 



Sk 



or 



with the normalisation constant Sk for each property i \ defined 



: Similar as in umbrella sampling, three different types of properties can be distirlbuished that 
< identity important states and should be included In the set 0 of properties. First, sfetes may be 
! important since they occur with large probability in the unbiased system {f{x) is luge). In this 
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c) of interest in 

I, there may 
in converged 
ily interested 
to ensure 
[d be added to 
irgence of the 



case, a property a(*)=l should be added to the set 0 of properties. Secon< 
fl$>portant since they contribute significantly to an average (a) of a property 
; the final analysis. Such a property should be added to the set © of properties, 
exist states that are not directly of interest but whose sampling is important to o| 

• estimates of other properties. For example/in molecular systems where one is 

• • iri lowenergy confonnations, sampling of high^nergy conformations is advantac 
t transitions between different regions of low energy. In this case, a function shoi 

the set © that is large for states that need to be sampled to increase com 

• estimates. Examples are given below. 

I. FittinaProcedurei 1: Fit to state with lame weighting factor 

; Procedure 1 fits a sampling distribution function such that it is particular large for Selected state 
with a large product fciCM- The description that follows assumes that|he sampling 
distribution function of Eq. (16) is used, that the centre ***** of the constraint isletermined by 
the fit, and that the value of the constant fc«~ is set to a positive value that is dttermined in.* 
series' of preliminary test runs. The procedure may be adapted, e.g., for cfcer 
j distribution functions. 

■ Based on the value of the product ft***), the state from the already sanfeled states Is 
• identified that has a comparatively large weighting factor. E g., the criterion to id^tify the state 
jto/ can be that it corresponds to the maximum 

i In Eq. (20), the maximum is taken over the previous simulations /. and ove|the states 

• i generated in each of the simulations; Nm=j-l is the number of simulations canid out so far. 
! The; centre Of the constraint is then set equal to the selected state, i.e., x,.«L==*w. The 
• seleetedstatex^ may also be used as the starting state for the simulation of the n|ct iteration. 



' The objective function that is maximised by procedure 1 is equal to 
. G= . max . max du . 
with the focal comparisons defined as 

and the states 36/ for the comparison set equal to xu . 
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Pitting Procedure' 2: Fit to states with small weighting factors. 



procedure 2 fits a sampling distribution function to a state with a particular mall weighting 
factor* The description assumes again that the sampling distribution function of ;q. (16) is used 
«Sijd that the centre xj.conxr of the constraint is determined by the fit The vafuefcf the constant 
Item** ; fs set to a negative value that may be determined in a series of preliminar test runs. The 
pfocedure may be adapted; e.g., for other sampling distribution functions with a defined centre. 
A. state xset is selected with a small product fatC/fas), e.g., the one that con [spondS to the 
minimum 

(23) 



TJj»e centre of the constraint xjjamar is set equal to the selected state x& . The ot^ctive function 
ttiat Is maximised by the procedure can be written as 

where the local comparisons are defined by Eq. (22). 



(24) 



More generally, a sampling distribution function can be fitted to a set of states 



Wth particular 



•small weighting factors, e.g., the centres xjj.c**r may be determined of the samp! «g distribution 



function of Eq. (18) with kjj,con*, set to negative values. For this purpose Nj. 
'■■ with 5=±\,2y. . JUicamr may be selected that correspond to the minima 

J^-V-^^l^^jpfe^ Jty^^^-^^y)^,))- (25) 
.. This* may be done by selecting first the state that corresponds to the minim 

with s^L , then select the state x**& with y=2 , and so on until s=Njf<mxr. The 
: ;of the constraint can then be set equal to the selected states, x**v . This procedi 
v theobjective ftinctlon of Eq. (24) with the local comparisons defined as 

.; Fitting Procedure 3: Linear least square fit. 



states x**te 




'Procedure 3 performs a linear least square optimisation to fit the vector vj of 
sampling distribution function to the product fh^Gfas] . The description 
sampling distribution function of Eq. (18) is used. The number of harmonic constra 



assun es 



para, ieters 



of the. 
that the 

tS, Nj,eotvtr f 
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ahd the centres, xjj^o- , of the harmonic constraints are determined before stAting the fit This 
can be done, for example, using a clustering algorithm. The strengths kjj,anatr m the constraints 
are obtained by a linear least square fit of the logarithm of the urvnorrridised sampling 



distribution function to the logarithm of the product pi^Oixi,^ of the weighting f 
the absolute value of the function o(x) . The objective function can be set equal 

with local comparisons defined by 

The procedure may be adapted, e.g., for other sampling distribution functions. 



ctors p.,, times 



(27) 



(28) 



Fitting Procedure 4: Fit to paira of states with large a nd small weighting factors, r&oectivelv. 

Procedure 4 fits the sampling distribution function to pairs of states x»n and j&fc with the two 
states being close together in the state space, and with p*h being relatively Jnall and p«n 
being relatively large. Such a pair can be identified, for example, using the criteril of Eq. (20) to 
select Xsm , and then select the state closest to It as . An alternative is t| select states 
generated in subsequent time steps with a Metropolis Monte Carlo Version 1, Ir a stochastic 
dynamics algorithm. The selected states are then used to define a bias away froji into the 
region around x»n , e.g., use the sampling distribution function of Eq. (16) with 
xjemotr = xjth +c(xw2 —jCtt/i ) , where c is a constant . 

The objective function maximised by this procedure is the same as the one liaximised by 
procedure 1, except that the state at which to evaluate the sampling distribution flnction should 
be defined as 3ft/=»,f+c(^-x»Ac/ 0 «), where *ac*«* denotes the state from the s|mp!ed states 
that is closest to xt,t . 

Sampling of conformations of molecules 

■ 

Sampling molecular conformations is on© of the main applications of moleci ar dynamics 
simulations. Examples are peptide and protein folding simulations, In which the c informational 
' space of a peptide or protein is investigated with the aim of identifying conforms ions that are 
stable under native conditions! or drug binding studies with the aim of idlntifying and 
characterising binding modes of drug candidates to their host proteins. As an ixampte, the 
following description describes how to apply the patent to the protein folding problen 
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Step 1. Perform a short molecular dynamics simulated annealing (S . Wrkpat 
Jr., M. P. Vecchi, "Optimization by Simulated Annealing", Sc/ence, 220, 4598, 
run starting from an extended conformation. The conformation, xomeei, res 
annealing run is used as the starting state of the simulation when executing 
time. In the first simulation, the system is restrained to the region around 
sampling distribution function of Eq. (15) is used with n,cmiu-=xamtai . The parar 
determines the strength of the harmonic constraint is set to a positive value 
based.on a set of preliminary test runs. 



bk, C. D. Gefatt 
[g83) 671-680.) 
itting from the 
2 for the first 
mat, i.e., the 
5ter kntotr that 
lat can be set 



starting state 



I (16). A single 



Step 2. Carry through a molecular dynamics simulation starting with the se 
(conformer), and the selected sampling distribution function pj{x) . 
Step 3;1. The weighting factors pv are determined with Eqs. (8-10). 
Step 3,2. Fitting Procedure 1 is used to fit the sampling distribution function of 
property is used to define the function 0(x) with Eq. (19) and to select a state wit* Eq. (20). The 
property is set equal to the inverse of the probability distribution of the energies If the system. 
The probability distribution of the energies Is estimated using Eq. (5) based on tie energies of 
the sampled conformations xv and their weighting factors p,/. The proper* is large for 
conformations with energies that are unlikely in the unbiased system, and 1 is small for 
conformations with energies that are likely in the unbiased system. Using this 
(ig) ensures that states are selected regardless of their energy, provided they 
that was relatively little sampled in previous simulations. Depending on the syst 
properties could be included in the set 0 of properties used to define the fiinctfen Oix). For 
example, for drug binding studies, the additional property might be set equal to thelwerse of the 
probability distribution of the position of the drug candidate relative to the prote|. This would 
ensure that different possible positions of the drug relative to the protein are sampl 
Step 4. Make 200 iterations of steps 2 to 4. 

Step 5. Determine the weighting factors pv. Use the weighting factors to estimate properties. 
e.g., the free energy as a function of the temperature, or the probability distributionlof the radius 
of gyration at a selected temperature. 



im 



arty in Eq. 
in a region 
i, additional 



Two-dimensional HP protein-folding model 



35 



This section applies the invention to study the two-dimensional HP protein-foldind 
compares it to simulated annealing and to multicanonical sampling. The conformatiJ 



model, and 
nal space is 
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Investigated of a chain consisting of 100 hydrophobic H and polar P monftnere. The chain 
moves on a two-dimensional grid. The chain is self-avoiding, i.e., no two monomers may occupy 
the same grid point ff two H monomers are In contact, they interact fevourablylgiving an energy 
. contribution of -2. There are no other interactions. The probability distributlorlof the system is 
5 set equal to /?(x)=texp(-£(x>2) where E(x) is the energy function. The problbffity distribution 
is not normalised. A Metropolis Monte Carlo algorithm is used to sample conlrmations of the* 
chain* The Metropolis Monte Carto algorithm uses a set of simple move* to change the 
conformation of the chain locally, one or two monomers per move. The samelset of moves is 
used In all the simulations, and no attempts were made to optimise the setpf Monte Carlo 
10 . moves. Fig. 5 shows a conformation. Hydrophobic H monomers that interactlfavourably if in 
contact are drawn as empty rectangles; polar P monomers are drawn as filled cimes. 

To study the conformational space using simple enumeration, some 17*1 0 4 ! conformations 
would have to be generated, which is beyond the capabilities of computers. IThe examples 
15 presented in this study change the conformation 3*10* times. The required calculations take a 
few hours on a personal computer. The properties that are assumed to be of merest and that 
are subsequently analysed are low energy states, states with small radii, the probability of states 
as- function of the energy and as a function of the radius. The first two properties Ilustrate global 
optimisations; the latter two properties illustrate state space integrals. I 

20 s | 

* The procedure used to study the folding of the grid model is similar to the one dAcribed above 
for the protein folding simulations. In particular, a sampling distribution functionlsimilar to the 
distribution function of Eq, (16) is used that applies a local constraint, andtthe sampling 
distribution function Is fitted using fitting procedure 1. In more detail: the constrainiig potential is 
25 defined based on the radius of gyration and on a comparison of the monomer coltacts present 
in.the current conformation x and the conformation selected by the fitting procesJx*/ . I.e., the 
. • sampling distribution function is set equal to 8 

pj(x}= fjCX^h(njbmed(x)±mrv^ . 129) 

In Eq. (29), n/orm^(x) Is equal to the number of monomer contacts formed i| the current 
30 .conformation xthat are not present in x*ei % fUnkai(x) is equal to the number fcf monomer 
contacts of x**i that are no more present in the current conformation x . R&r(x) islhe radius of 
gyration of conformation x , and R*/ is set equal to the radius of gyration of h*t . Since a 
Metropolis Monte Carto algorithm (version 1) is used to generate the conforthations, the 
normalisation constant fj in Eq. (29) needs not to be known. The constants h ar|J are set 
35 equal to 0.003 and 0.5, respectively. The constants were determined in a series of t|st runs. 
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tpoF^he simulations (step 2), a simple Metropolis Monte Carlo algorithm (vers* 
■Monte Carlo algorithm generates new conformations by randomly changing! 
^fngfe monomer or two adjacent monomers at a time. New conformations I 
.tejested using the criterion of Eq. (6). A single run of the simulation algorithm I 
.of 3'000'QOO Monte Carlo moves. The first 500*000 moves are consider 
phase, and are not used Then, 50 conformations are taken every 50*000 mov 
50 cehformatfons xjj for each of the simulations / 100 simulations am carried 
of 3*10 8 states are generated, of which 5'GOO are used for the analysts, e.jgr M tc 
cfi?tributiorTfunctionslnstep3. . 



1) is used. The 
le position of a 
accepted or 
[step 2, consists, 
equilibration- 
This results in 
Thus, a total 
itthe.sampling 
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Two iaroperties are used to define the function 0(x) with Eq. (19) and to 
dfetribution function using procedure 1. The first property is the inverse of I 
distribution of the energies of the system. The probability distribution of the energl 
using. Eq. (5) based on the energies of . the sampled conformations *v and 
factor^ fa . The property is large for conformations with energies that are 
unbiased system described by p{x), and it Is small for conformations with enJ 
likely in the unbiased system. The property ensures that states are selected regj 
eneigy*tf they are in a region that was little sampled in previous simulations. 

The second property used in Eq. (19) is the fnveree of the probability distribution 
gyjSfion of the system. This property ensures that states are selected regardless 
. ofgyrason, if they are in a region that was little sampled in previous simulations, 
: the, properties to use for the fitting procedure is compatible with the set of prop 
assumed to be of interest (low energy states, states with small radii, the probabifil 
. function of the energy and as a function of the radius; see above). 
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Results from five independent test runs are analysed. Fig. 6 plots the probability <| states as a 
• funetion of the energy. The agreement between the results of the different runs isfcood. Some 
^differences exist forthe estimates of the probabilities of states with low energies. Thl differences 
jn the estimates suggest that there exist different kinds, of low energy states, and |iat the runs 
diftet in the set of low energy states sampled. Low energy conformations with enerdfes around - 
iSO'ere rare, i.e., they are estimated to occur with a probability of about 10" 46 . 

Vig. ? illustrates the cdnvergence of the estimates of the probability of low energjlstates from 
Fig; 6. Shown are the estimates of the probability as a function of the progress of thl simulation. 
.Low, energy conformations are rare, and finding them is difficult. They are found onrlatter some 
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initial sampling that explores the state space. Once, a first state with given Inergy has been 
found, the estimate of the probability converges rapidly, and does remain stabl| over the rest of 
the run. 

Fig. 8 plots the energy of the states from Run B as a function of the progress oflhe run. The run 
startsby sampling high energy states. As the run progresses, it finds states of Increasingly low 
energy. The system gets never trapped After having found a state of low Inergy. the run 
contiriues by sampling new regions of the state space, and eventually finds a |ew low energy 
state. 

.Fig. 9 shows estimates of the radius distribution of states. The estimates from tit different runs 
are all smooth functions and agree quantitatively. 

Comparison with the sta te of the art 

In summary, the five runs are successful in identifying low energy and compact sj 
quantitative estimates of the probability of the energy of states and of the radii of 
test runs are compared to five simulated annealing runs and to five runs of 
sampling. The same Metropolis Monte Carlo algorithm is used in all the runs, ai 
runs, the conformations are updated 3*10 8 times. In the simulated annealing run ( 
C. D. Gelatt Jr., M. P. Vecchi, "Optimization by Simulated Annealing", Sciem 
(1983) 671-680.). the temperature is lowered linearly in 5000 steps from 1.5 
multicanonical runs (Bemd A Berg and Thomas Neuhaus "Multicanonical ens 
•approach to simulate first-order phase transitions", Phys. Rev. Let.. 68 (1992| 
umbrella potential is updated 100 times. 



:, and give 
tes. The five 
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as in the test 
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Table 1 compares the runs with respect to their performance for the optimisation Iroblem "find 
low energy states". The lowest energy found in different runs varies between -fo and -94. 
From all the runs, the test runs with the present invention identify the best state (en|gy of -94). 
.In average, the test runs find states with an energy of -88. This compares favojrably to the 
average of -66.4 from the multicanonical runs, and to the average of -63.6 from tie simulated 
annealing runs. 



Min Score 



Adumb 

Annealing 

Test 



Run A Run B Run C 



-86 
-86 
-86 



-86 
-88 
-90 



-86 
-84 
-90 



Run D 

"3>6 



Run E 



-80 
-80 



I -88 
I -80 
-94 
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Table 1 : Lowest energies found in each of the runs. Data is shown for the five ast 
runs (Test)', the five multicanonical or adaptive umbrella sampling runs (Adum ), 
and the five simulated annealing runs (Annealing). 



Runs with the present invention do not get trapped (see Fig. 8). This is different 
situation with simulated annealing and mutticanonical simulations. In simulated 
(see Fig. 10), the simulations start at a high temperature at which transitions bet 
regions occur frequently. As the temperature Is lowered, sampling is increasingly 
region around a low energy state. If different regions with comparable low energl 
the present system - several simulated annealing runs are required to identif 
regions. 



Multicanonical simulations of high dimensional system have the tendency to 
11) close to the first low energy state found. Although multicanonical simulations 
gu&rantee that low and high energy states are sampled with comparable probabilif 
guarantee that transitions between the different regions occur frequently. 



Fig. 12 shows estimates of the probability of states with given energy from the 
runs and compares them to the corresponding estimates from the test runs 
average of the estimates from the test runs Is shown as continues line in Fig. 12. 
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nuIHcanonical 
(Fig. 6). The 
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from all the runs agree In general. The estimates for high-energy states agre > well. Some 
differences, are evident at low energies. At low energies, the probabilities estimate from the test 
runs are larger than those from the multicanonical runs. This indicates that the test runs find 
more low energy states. 1 

The simulated annealing runs are irrelevant in such a comparison of equilibriun properties of 
the system, since no method has so far been described to derive converge! state space 
integrals from simulated annealing runs. 

Value at Risk of a Financial Portfolio 

The value at risk of a financial portfolio gives an estimate of the future value lof a financial 
portfolio in adverse market situations. It is equal to a percentile of the probabilityldistribution of 



the future value of the portfolio. E.g. t the 99% value at risk is equal to the futui 



portfolio that separates the worst 1% from the best 99% of the possible future values of the 



portfolio. Efficient estimation of the value at risk is challenging since it requires sar 



value of the 
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rare state, i.e., particular adverse market situations, and the location of the rah 
general not known before starting the simulation? and must be found. 



states is. in 



Step 1. Select the distribution function of the (unbiased) system as the in ial sampling 
distribution function, i.e., pi{x)=p(x). A Monte Carlo algorithm is used for the sir* illations. The 
Monte Carlo algorithm produces states independently of any starting state, Thi s, no starting 
state is needed. A Monte Carlo algorithm must exist to produce un-correlated sta is distributed 
according to pi{x). This restricts the set of sampling distribution functions that can e used. 
Step 2. Carry through a Monte Carlo simulation using the selected sampling distrit iHon function 

Step 3.1. The weighting factors pu can be determined with Eqs. (8) and (10). The|iormalisation 
•constants ff J are known and need not to be determined with Eq. (9) 



Step 3.2. Fitting Procedure 1 is used to fit the sampling distribution tundion of Eq 
property Is used to define the function C(x) with Eq. (19) and to select a state with 
property Is defined as follows. First, determine a limit identifying large losses, e.j 
occur with a probability of less than 10%. Then, set the property equal to a functiqji that is 0 for 
small losses and 1 for large losses. 



15). A single. 
=q. (20). The 
, losses that' 



Step 4. Make 200 iterations of steps 2 to 4. 
Step 5. Estimate the weighting factors pi,. Use the weighting factors to djtermine the 
probability distribution of the value of the portfolio, and determine the value at risk. 
An alternative would be to use the sampling distribution function defined by |q. (18), the 
Metropolis Monte Carlo Version 2 algorithm, and the fitting procedures 2 or 3. 
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Claims 



sa 
a- fitting 



A Method for sampling a state space by iteratively generating states and heir 

weighting factors py , wherein the Index / Is the iteration parameter and the idex t 

distinguishes different states generated by an iteration /, the method Is a isisting of 

* first step for selecting ah initial sampling distribution function pi(x), 

a fifth step for performing an analysis 
and an iteration procedure including 

a second step for generating Nj states xjj by a numerical sampling algorith nand 
a fourth step for testing at least one criterion to decide whether to continue th > iteration 
procedure or to stop the iteration procedure and to go to a fifth step in order t . perform the 
analysis using the simulated data, 

characterized in that the iteration procedure further includes 
a third step determining weighting factors p., for states xt/ generated so far >y using 
tmpling distribution functions p f (x) determined so far and 
ting step for determining a sampling distribution function p s (x) for the ne t iteration by 
fitting Pj(x) to ptjOfaj) for states xij generated so far, wherein Ofas) iss function, 
respectively a property, of the states Xijt . 

Method as claimed in daim 1, wherein the sampling distribution function pj{x i of at least 
one iteration is fitted such that it maximises an objective function preferably d< fined as a 
function of local comparisons between the sampling distribution function and fie product 

A/0(*/,r). 

Method as claimed in claim 1 or 2, wherein the sampling distribution function p/( x ) of at 
least one iteration is fitted such that the sampling distribution function pj{x) i Jlarge for at 
least one state *v with a large product pvO(xi,<) , and tends to be small for sptes with a 
small product pvG{pS)- 

Method as claimed in one of claims 1 to 3, wherein the sampling distribution f notion 
pj(x) of at least one iteration is a function with at least one constraint, prefer? )ly having 
an extreme value in the region of at least one selected state xj^^r and havir g vanishing 
values at a distance from the at least one selected state xj wC o>wr , wherein the < onstraint is 
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7. 



preferably a harmonic constraint with at least one constant kcow defining th| strength of 
the constraint at the selected state xj,c 



Method as claimed In claim 4, wherein the sampling distribution function pj{ :) of at least 
one iteration is the distribution function of the system with, constraints or a mpticanonical 
distribution function with constraints. 

Method as claimed in claim 1 or 2, wherein the numerical sampling algorithm of at least 
one iteration generates correlated states xp , wherein the sampling distribut m function 
pj(x) preferably has a maximum that biases sampling into a sub-region oft e state 
space, wherein the sampling distribution function pj(x) is preferably fitted sjch that the 

maximum is in a region where the product p^Ofaj) has a maximum, and u lerein the 
sampling preferably starts from a state close to the maximum of the samplinf distribution 
ftjnction pj(x). 



don 



Method as claimed in one of claims 4 to 6, wherein the fitting of pj(x) is 
states x u for which the product j5uO(x*.r) has extreme values and by using 
Xij to define the region in which pj(x) has extreme values. 



by selecting 
he selected 



10 8. • Method as claimed in one of claims 1 to 6, wherein parameters of the sampfi tg distribution 
function pj(x) of at least one iteration are determined by a linear least square fit of the 
logarithm of the unformalized sampling distribution function pj(x) to the lo^rithm of the 
product pjOipS) 

15 9. Method as claimed in one of claims 1 to 8, wherein the normalization constalt of the 

sampling distribution function p/x) of at (east one iteration is estimated fronflthe sampled 
' ' states x M and their weighting factors p?/ . . 

10. Method as claimed in one of claims 1 to 9. wherein at least three iterations ai& done. 



1 1 . Method as claimed in one of claims 1 to 10, wherein the function d(x) is a ft ictiori of a 
set of at least two functions ©=^a(x)Qz(x). . „0jw(x)} , preferably including |t least one 
of the following functions 
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at least one property for which at least one estimate is derived in the analysis of the fifth 
step, 

at least one function that Is larae for states that must be sampled to ensure ti msitions 
between important regions, 
5 ' the inverse of the probability distribution function of at least one property of tl j system, 

and 

the inverse of the probability distribution of the negative logarithm of the distri lution 
function of the system . 

o 12. A computer software product for sampling a state space by iteratively genera ing states 
and their weighting factors /5».r , wherein the Index /is the iteration param terandthe 

index X distinguishes different states *m generated by an iteration /, said predict is 
characterized by a computer-readable medium in which program instructionsfire stored, 
which instructions, when read by a computer, enable the computer to 
s select an Initial sampling distribution function p\(x) , 

to execute a fifth step for performing an analysts 
and to execute an iteration procedure including 

a second step for generating Nj states x» by a numerical sampling algoritrAi and 

a fourth step fortesting at least one criterion to decide whether to continue th( iteration 
0 procedure or to stop the iteration procedure and to go to a fifth step in order t( perform the. 

analysis using the simulated data, 
• characterized in that the Iteration further includes 

a third step determining weighting factors pv for states xt, generated so far »y using 

sampling distribution functions p,{x) determined so far and 
s a fitting step for determining a sampling distribution function pj(x) for the ne> iteration by 

fitting p/x) to pvOfa) for states xu, generated so far, wherein 0{x») is a|inction, 

respectively a property, of the states xv ■ 
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The present Invention is an incremental umbrella sampling method to improve thl performance 
of established sampling methods. It is sampling the state space by iteratively gei brating states 
xv and their weighting factors pu by fitting the sampling distribution function pA c) of the next 
iteration to at least one weighted property of the already sampled states. This m« ns that pj(x) 
is fitted to the product j&Ofci), in which the p.., are the weighting factore ar I G{xi,) is a 
function respectively a property of the states jg, . The number of states xv and te number of 
weighting factors pv Is incremented with each iteration. In order to have a coi sistent set of 
weighting factors p fJ , the weighting factors are recalculated in each iteration for a • respectively 
for a set of selected, states. By fitting pj(x) in the state spade it is possible 3 use all the 
Information of pv and Cfa) for the states x„ generated so far. The fitting step allows to use 
different fitting strategies. For example the fitting can bias the sampling away fron areas where 
Intensive sampling has been done in the preceding iterations, or the sampling cs i be directed 
along local gradients respectively towards local minima or maxima of one or aev jral weighted 
properties. In each of the iterations, the sampling distribution function is fitted in a way to 
improve the overall sampling of the state space. The method supports I lUlthobjective 
optimisations. State space integrals can be solved. It reduces the probability that he system is 
trapped. The invention is general. It can be used with different sampling method? : in particular 
with Monte Carlo sampling, Metropolis Monte Carlo sampling, or dynamic simulati ns. It can be 
combined with the concepts of simulated annealing and muKicanonica! sampling, it provides a 
general framework that can be adapted to the system and the observables of interest 
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